Q1_projects_under_7days

Question #1

For burn projects that were wrapped up in less than a week:

a.) How well did the number of days for which thermal anomalies were detected corroborate with the number of days that burners indicated in their project reporting?

b.) Is there an apparent detection limit (e.g., acres burned or some other metric) under which the satellite does not seem to perform well in terms of picking up thermal anomaly? How does the finding different for broadcast burn vs. pile burn projects? How does the apparent detection limit compare to the US EPA’s Rx burn reporting thresholds (50-acre for broadcast burn; 25-acre for pile burn)?

Approach

• Using CalMAPPER data, keep only treatments that 1) included Rx burns, 2) lasted 7 or fewer days
• Filter out all VIIRS points outside CalMAPPER polygons
• for each polygon, sum together the number of unique days with CalMAPPER activity
• Do the same for VIIRS points within each polygon.
• Calculate the difference in duration and look at the distribution.

Import data

# VIIRS fire detects
viirs <- read_rds(file.path(ref_path, 'VIIRS/CA_2020/viirs_nonAg.rds')) 

# CalMAPPER data
CM <- read_rds(file.path(ref_path, 'CalMAPPER/treatment-activities.rds'))

Prep CalMAPPER data

Just get burn data out of CalMAPPER.

# check if cultural burning can be either pile burn or broadcast. all broadcast this year at least.
CM %>% 
  as_tibble %>% 
  filter(ACTIVITY_DESCRIPTION == 'Cultural Burning') %>% 
  distinct(ACTIVITY_DESCRIPTION, ACTIVITY_NAME)

## # A tibble: 2 × 2
##   ACTIVITY_DESCRIPTION ACTIVITY_NAME                             
##   <chr>                <chr>                                     
## 1 Cultural Burning     Prescribed Fire - Cultural Burning        
## 2 Cultural Burning     Prescribed Fire - Oak Woodland Restoration

# CalMAPPER -- <= 7 days, trts include burns
CM_burn <- CM %>% 
  # create columns for burn activity
   mutate(burn = ACTIVITY_DESCRIPTION %in% 
           c('Broadcast Burn', 'Cultural Burning', 'Pile Burning'),
          pile_burn = ACTIVITY_DESCRIPTION == 'Pile Burning', 
          broadcast_burn = ACTIVITY_DESCRIPTION %in% 
            c('Broadcast Burn', 'Cultural Burning'),
         .after = AQ_ID) %>% 
  # get duration of burn activities. +1 to duration (e.g. 2/1-2/1 = 0. It should be 1 day.) 
  mutate(duration = difftime(ACTIVITY_END, ACTIVITY_START, units = 'days') + 1) %>% 
  # include only burn activity
  filter(burn == T) %>% 
  # get area of treatment too. might be useful
  mutate(trt_acres = units::set_units(st_area(.), 'acres') ) %>% 
  # keep only the useful columns
  select(TREATMENT_ID, AQ_ID, burn:broadcast_burn, TREATMENT_NAME, ACTIVITY_NAME, ACTIVITY_DESCRIPTION, ACTIVITY_START, ACTIVITY_END, duration, ACTIVITY_STATUS, TREATED_ACRES, trt_acres, BROAD_VEG_TYPE, AGENCY_NAME, FY)

Just look at projects from 2020. CalMAPPER is supposed to be complete only FY19/20 and up. I downloaded VIIRS data from Jan1-Dec31. Check out the projects that span calendar years in and out of 2020… If project spans calendar year, the start date is always on Dec 31st. Seems like a default value.

Update: Dec 31 1600 PST is the same as Jan 1 0000 UTC. It seems like a default value whenever the start date isn’t known. So, remove the 2021 start dates, but 12/31/2019 starts are ok.

# projects spanning 2020 calendar year
CM_burn %>% 
  filter(year(ACTIVITY_END) > 2020 & year(ACTIVITY_START) == 2020 |
           year(ACTIVITY_START) < 2020 & year(ACTIVITY_END) == 2020) %>% 
  pull(ACTIVITY_START) %>% 
  unique

## [1] "2019-12-31 16:00:00 PST" "2020-12-31 16:00:00 PST"

# filter just 2020
CM_burn2020 <- CM_burn %>% 
  filter(year(ACTIVITY_END) == 2020) 

Let’s also just look at the distribution of start and end dates. Are their common dates that the records tend to show? Spikes in the figure would indicate inaccurate dates.

CalMapper_dates <- cowplot::plot_grid(
  ggplot(CM_burn, aes(ACTIVITY_START)) +
  geom_bar(),
ggplot(CM_burn, aes(ACTIVITY_END)) +
  geom_bar(),
nrow = 2
)
CalMapper_dates

ggsave(plot = CalMapper_dates, 
       filename = 'tmp_outputs/CalMAPPER_activity_dates.png', width = 6, height = 4)

This data set has 2 treatments that have the exact same geometries, even though they have different ID numbers. This will cause issues when trying to count the number of unique fire detect days within the polygons. e.g. one treatment has a burn in january and the other treatment has burns from feb-may. The additional detects will seem anomolous for the first treatment, when in fact it should be part of the second treatment.

Since there are only 2 treatments that are duplicated, manually merge them. In the future, manual cleaning up of the data should be avoided.

# calmapper treatment ids. 
CM_trtIDs <- CM_burn2020 %>% 
  group_by(TREATMENT_ID) %>% 
  summarise()

# find which rows have duplicates
mat_eq <- st_equals(CM_trtIDs, sparse = F, remove_self = T)
duplicated_trts <- CM_trtIDs$TREATMENT_ID[apply(mat_eq, 1, sum) == 1] 
CM_burn2020 %>% filter(TREATMENT_ID %in% c(duplicated_trts))

## Simple feature collection with 22 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13463530 ymin: 4732515 xmax: -13461520 ymax: 4737232
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
##    TREATMENT_ID AQ_ID burn pile_burn broadcast_burn            TREATMENT_NAME
## 1          9553 60497 TRUE      TRUE          FALSE           N-05-19 Phase I
## 2         14544 62420 TRUE      TRUE          FALSE North Fork American River
## 3         14544 63060 TRUE      TRUE          FALSE North Fork American River
## 4         14544 67008 TRUE      TRUE          FALSE North Fork American River
## 5         14544 68189 TRUE      TRUE          FALSE North Fork American River
## 6         14544 74051 TRUE      TRUE          FALSE North Fork American River
## 7         14544 74066 TRUE      TRUE          FALSE North Fork American River
## 8         14544 74898 TRUE      TRUE          FALSE North Fork American River
## 9         14544 77458 TRUE      TRUE          FALSE North Fork American River
## 10        14544 78102 TRUE      TRUE          FALSE North Fork American River
##                         ACTIVITY_NAME ACTIVITY_DESCRIPTION      ACTIVITY_START
## 1  2019 CNG & Fuels Crew Pile Burning         Pile Burning 2020-01-05 16:00:00
## 2               2020 CNG Pile Burning         Pile Burning 2019-12-31 16:00:00
## 3               2020 CNG Pile Burning         Pile Burning 2020-01-19 16:00:00
## 4               2020 CNG Pile Burning         Pile Burning 2020-01-26 16:00:00
## 5               2020 CNG Pile Burning         Pile Burning 2020-02-02 16:00:00
## 6               2020 CNG Pile Burning         Pile Burning 2020-02-09 16:00:00
## 7               2020 CNG Pile Burning         Pile Burning 2020-02-23 16:00:00
## 8               2020 CNG Pile Burning         Pile Burning 2020-03-01 16:00:00
## 9               2020 CNG Pile Burning         Pile Burning 2020-03-08 17:00:00
## 10              2020 CNG Pile Burning         Pile Burning 2020-03-15 17:00:00
##           ACTIVITY_END duration ACTIVITY_STATUS TREATED_ACRES        trt_acres
## 1  2020-01-09 16:00:00   5 days        Complete           1.5 1501.376 [acres]
## 2  2020-01-16 16:00:00  17 days        Complete           6.5 1501.376 [acres]
## 3  2020-01-23 16:00:00   5 days        Complete           7.0 1501.376 [acres]
## 4  2020-01-30 16:00:00   5 days        Complete           6.5 1501.376 [acres]
## 5  2020-02-06 16:00:00   5 days        Complete           6.0 1501.376 [acres]
## 6  2020-02-20 16:00:00  12 days        Complete           8.0 1501.376 [acres]
## 7  2020-02-27 16:00:00   5 days        Complete           3.0 1501.376 [acres]
## 8  2020-03-05 16:00:00   5 days        Complete           4.0 1501.376 [acres]
## 9  2020-03-12 17:00:00   5 days        Complete           4.0 1501.376 [acres]
## 10 2020-03-19 17:00:00   5 days        Complete           1.0 1501.376 [acres]
##    BROAD_VEG_TYPE AGENCY_NAME    FY                          SHAPE
## 1          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 2          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 3          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 4          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 5          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 6          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 7          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 8          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 9          Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...
## 10         Timber    CAL FIRE 19/20 MULTIPOLYGON (((-13461946 4...

# change the duplicated treatment ID
CM_burn2020$TREATMENT_ID[CM_burn2020$TREATMENT_ID == 9553] <- 14544

# rerun to update
CM_trtIDs <- CM_burn2020 %>% 
  group_by(TREATMENT_ID) %>% 
  summarise()

Export this data so that you can use it for the other questions.

write_rds(CM_burn2020, file.path(ref_path, 'CalMAPPER/burnactivities_2020.rds'))
write_rds(CM_trtIDs, file.path(ref_path, 'CalMAPPER/treatmentgeometries_2020.rds'))

Since I’m going to filter the VIIRS detections based on the CalMAPPER polygons, I should create a dataset with polygons with an extra buffer, where the radius is half that of the VIIRS resolution. This is just in case there’s a VIIRS detection right outside the polygon–I don’t want to miss that. And export.

# buffer polygons
st_crs(CM_burn2020)$units # make sure buffer will be in meters

## [1] "m"

CM_burn2020_buffered <- st_buffer(CM_burn2020, 375/2) 
CM_trtIDs_buffered <- CM_burn2020_buffered %>% 
  group_by(TREATMENT_ID) %>% 
  summarise()
# export
write_rds(CM_burn2020_buffered, file.path(ref_path, 'CalMAPPER/burnactivities_2020_buffered.rds'))
write_rds(CM_trtIDs_buffered, file.path(ref_path, 'CalMAPPER/treatmentgeometries_2020_buffered.rds'))

Finally, keep treatments where burn activity was <= 7 days.

CM_7days <- CM_burn2020_buffered %>% 
  filter(duration <= 7)

Prep VIIRS data

Filter out VIIRS data outside of the the CalMAPPER polygons. I’m using the buffered CalMAPPER polygons.

# get the viirs points contained in the treatment polygons
viirs_t <- st_transform(viirs, st_crs(CM_7days))

# calmapper treatment ids. 
CM_7days_trtIDs <- filter(CM_trtIDs_buffered, TREATMENT_ID %in% unique(CM_7days$TREATMENT_ID))
# CM_7days_trtIDs_unbuff <- filter(CM_trtIDs, TREATMENT_ID %in% unique(CM_7days$TREATMENT_ID))
# mapview(CM_7days_trtIDs) + mapview(CM_7days_trtIDs_unbuff)


# find those that intersect
viirs_int <- viirs_t %>% 
  filter(CONFIDENCE != 'l') %>% 
  st_intersection(CM_burn2020_buffered) %>% 
  mutate(rowname = rownames(.), 
         pointID = floor(as.numeric(rowname)),
         overlap_trts = grepl('\\.', rowname),
         .before = LATITUDE) %>% 
  group_by(pointID) %>% 
  mutate(overlap_trts = any(overlap_trts == T)) %>% 
  ungroup()

Find if any points sit on overlapping treatment. It looks like a handful of VIIRS points intersect overlapping polygons. This seems to often happen when there’s burn prep ahead of a broadcast burn. Out of the 1463 VIIRS points that intersect the CalMAPPER polygons, only 101 (7% of detects) are over overlapping treatments. Ignore for now and add a caveat for potential artifacts?

# find if any points sit on overlapping treatment
viirs_overlap <- filter(viirs_int, overlap_trts == T)

# number of viirs points vs those that overlap polygons
viirs_int %>% 
  as_tibble() %>% 
  distinct(pointID, .keep_all = T) %>% 
  count(overlap_trts)

## # A tibble: 2 × 2
##   overlap_trts     n
##   <lgl>        <int>
## 1 FALSE          721
## 2 TRUE          1404

# visualize
mapview::mapview(viirs_overlap, cex = 3, zcol = 'TREATMENT_ID', legend = F) + 
  mapview::mapview(CM_7days, zcol = 'TREATMENT_NAME', legend= F, alpha.regions = 0.3)

Q1a

How well did the number of days for which thermal anomalies were detected corroborate with the number of days that burners indicated in their project reporting?

For every treatment polygon, calculate the total duration of burning activity and number of unique days for intersecting VIIRS detections. Compare the two calculations.

I also want to take a quick look at timing. For each treatment polygon, I calculated the mean julian day of both the VIIRS and CalMAPPER dates and compared the two. Gives you a general idea of temporal mismatch. For a more careful analysis, I might estimate the probability of a VIIRS detection in the treatment within 2 days of the activity. But can do that at a later point.

# duration + timing of CalMAPPER
CM_duration_timing <- CM_7days %>% as_tibble %>% 
  mutate(mean_jday = (yday(ACTIVITY_START) + yday(ACTIVITY_END))/2) %>% 
  group_by(TREATMENT_ID) %>% 
  summarise(CM_days = as.integer(sum(duration)), 
            CM_mean_jday = mean(mean_jday))

CM_vs_viirs_7d <- viirs_int %>% as_tibble() %>% 
  group_by(TREATMENT_ID) %>% 
  summarise(n_viirs = n(),
            viirs_days = length(unique(jday)), 
            viirs_mean_jday = mean(jday)
            ) %>% 
  right_join(CM_duration_timing) %>% 
  mutate(days_diff = CM_days - viirs_days,
         days_diff_abs = abs(days_diff), 
         jdays_diff = CM_mean_jday - viirs_mean_jday,
         jdays_diff_abs = abs(jdays_diff))

Visualize: 1) How treatment duration and timing differs between CalMAPPER records and VIIRS detections.

# distribution of viirs vs CM # of burn days. seems like in general little bias, with some treatments with way more activity recorded in CalMAPPER
p1 <- ggplot(CM_vs_viirs_7d, aes(days_diff)) +
  geom_rect(aes(xmin = -Inf, xmax = 0, ymin = -Inf, ymax = Inf), fill = "lightcyan2") +
  geom_rect(aes(xmin = 0, xmax = Inf, ymin = -Inf, ymax = Inf), fill = "indianred1") +
  geom_histogram(color = 'grey20', fill = 'grey70', alpha = .9) +
  geom_vline(xintercept = 0, lty = 2) +
  labs(x = "Difference between treatment duration \n(CalMAPPER - VIIRS)", 
       y = "# of treatments") +
  scale_x_continuous(n.breaks = 10)

# look at the absolute differences
q90 <- quantile(CM_vs_viirs_7d$days_diff_abs, .9, na.rm = T) 
p2 <- ggplot(CM_vs_viirs_7d, aes(days_diff_abs)) +
  geom_histogram(color = 'grey20', fill = 'grey70', alpha = .9, binwidth = 1)+
  geom_vline(xintercept = q90, lty = 2) +
  geom_text(aes(label = glue::glue("90th percentile of treatments \n(i.e. Between datasets, discrepency in duration \nusually under {q90} days)")), 
            x = q90 + 1, y = 15, size = 2.5, hjust = 0) +
  scale_x_continuous(n.breaks = 10) +
  labs(x = 'Absolute difference in treatment duration \n(|CalMAPPER - VIIRS|)',
       y = '# of treatments')

cowplot::plot_grid(p1, p2, rows = 1, scale = .85) +
  cowplot::draw_text('Number of burning days: CalMAPPER vs. VIIRS', y = .95)

Duration might look ok, but what about timing? Not great! And these are the treatments with detections too–plenty without any VIIRS points showing up.

# proportion of treatments with a difference of 2 days or less
mean(CM_vs_viirs_7d$jdays_diff_abs <= 2, na.rm = T)

## [1] 0.2962963

# differences in mean day of the year
p1 <- ggplot(CM_vs_viirs_7d, aes(jdays_diff)) +
  geom_rect(aes(xmin = -Inf, xmax = 0, ymin = -Inf, ymax = Inf), fill = "lightcyan2") +
  geom_rect(aes(xmin = 0, xmax = Inf, ymin = -Inf, ymax = Inf), fill = "indianred1") +
  geom_histogram(color = 'grey20', fill = 'grey70', alpha = .9) +
  geom_vline(xintercept = 0, lty = 2) +
  labs(x = "DDifference between treatment's avg julian day \n(CalMAPPER - VIIRS)", 
       y = "# of treatments") +
  scale_x_continuous(n.breaks = 10)

# look at the absolute differences
q90 <- quantile(CM_vs_viirs_7d$jdays_diff_abs, .9, na.rm = T) 
q50 <- quantile(CM_vs_viirs_7d$jdays_diff_abs, .5, na.rm = T) 
p2 <- ggplot(CM_vs_viirs_7d, aes(jdays_diff_abs)) +
  geom_histogram(color = 'grey20', fill = 'grey70', alpha = .9, binwidth = 1)+
  geom_vline(xintercept = c(q50, q90), lty = 2) +
  geom_text(aes(label = glue::glue("90th percentile \n(i.e. ~{round(q90)} days)")), 
            x = q90 + 5, y = 15, size = 2.5, hjust = 0) +
  geom_text(aes(label = glue::glue("50th percentile \n(i.e. ~{round(q50)} days)")), 
            x = q50 + 5, y = 15, size = 2.5, hjust = 0) +
  scale_x_continuous(n.breaks = 10) +
  labs(x = 'Absolute difference in treatment duration \n(|CalMAPPER - VIIRS|)',
       y = '# of treatments')

cowplot::plot_grid(p1, p2, rows = 1, scale = .85) +
  cowplot::draw_text('Avg julian day for burns: CalMAPPER vs. VIIRS', y = .95)

Visualize: 2) how many polygons are without any VIIRS detections. 60%!

nodetects <- mean(is.na(CM_vs_viirs_7d$viirs_days))
nodetects

## [1] 0.604878

Visualize: What about the polygons with no detections–what’s the distribution of detections within treatment? Are there any patterns associated with these polygons–Could size of treatment matter? Burn type?

# add a column for VIIRS detections
CM_7days2 <- CM_7days %>% 
  left_join(select(CM_vs_viirs_7d, TREATMENT_ID, n_viirs)) %>% 
  mutate(VIIRS_detects = !is.na(n_viirs),
         n_viirs = replace_na(n_viirs, 0), .after = n_viirs)

# look at number of detections within trts?
CM_7days2 %>% 
  as_tibble %>% 
  distinct(TREATMENT_ID, .keep_all = T) %>% 
  ggplot(aes(n_viirs)) +
  geom_histogram(color = 'grey20', fill = 'grey70', alpha = .9, binwidth = 1)

# prep data
CM_size_burntype <- CM_7days2 %>% 
  as_tibble() %>% 
  group_by(TREATMENT_ID) %>%
  mutate(any_broadcast_burn = any(broadcast_burn == T)) %>%
  distinct(TREATMENT_ID, trt_acres, VIIRS_detects, any_broadcast_burn) %>% 
  ungroup

Look at burn type’s influence on detections–slightly more in treatments with broadcast burns than pile burns. But the majority of either treatment are not associated with VIIRS detections.

# look at burn type vs prob of detection
count(CM_size_burntype, any_broadcast_burn, VIIRS_detects) %>% 
  mutate(any_broadcast_burn = ifelse(any_broadcast_burn == T, 'Broadcast burn', "Pile Burn"), 
         VIIRS_detects = ifelse(VIIRS_detects == T, 'Detected', "Not detected")) %>% 
  ggplot(aes(any_broadcast_burn, n, group = VIIRS_detects)) +
  geom_col(aes(fill = VIIRS_detects), 
           position = position_dodge2()) +
  labs(x = 'Treatment type', 
      fill = "VIIRS", 
      y = "# of treatment polygons") +
  scale_fill_brewer(palette = 'Set1', direction = -1) +
  theme(legend.position = 'bottom')

# slightly more detections with broadcast burns than pile burns
CM_size_burntype %>% 
  group_by(any_broadcast_burn) %>% 
  summarise(perc_detections = mean(VIIRS_detects))

## # A tibble: 2 × 2
##   any_broadcast_burn perc_detections
##   <lgl>                        <dbl>
## 1 FALSE                        0.333
## 2 TRUE                         0.468

Ok, well is that bc broadcast burns might be over a larger area than pile burns? Not necessarily. If anything, the importance of treatment size seems less for broadcast burns. This figure nonetheless shows that VIIRS generally wont do a good job of detecting activity until >1000 acres, much greater than the thresholds of 25 and 50 acres for burn piles and broadcast burns (dashed vertical lines), respectively. Meanwhile, the median treatment is 97 acres.

# look at size of treatment on probability of any VIIRS detections. 
CM_size_burntype %>% 
  mutate(burntype = ifelse(any_broadcast_burn == T, 'Broadcast burn', 'Pile burn')) %>% 
  ggplot( aes(as.numeric(trt_acres), as.numeric(VIIRS_detects), 
                             group = burntype, color = burntype)) +
  geom_jitter(width = 0, height = .00, alpha = .3) +
  geom_vline(data = data.frame(burntype = c('Broadcast burn', 'Pile burn'), x = c(50, 25)), aes(xintercept = x, color =burntype), lwd = 1, lty = 2) +
  geom_smooth(method="glm", method.args = list(family = "binomial")) +
  labs(x = 'Treatment size (acres)', y = "Probability of any VIIRS detects",
       title = 'VIIRS detections vs. size', 
       subtitle = "*EPA's reporting thresholds shown by vertical lines", 
       color = 'Burn type') +
  #scale_x_continuous(n.breaks = 10) +
  scale_x_log10(n.breaks = 10) +
  scale_y_continuous(n.breaks = 10) +
  scale_color_brewer(palette = 'Set1', direction = -1) +
  theme(legend.position = 'bottom')

# whats the typical size of treatments? 
median(CM_size_burntype$trt_acres)

## 97.33222 [acres]

ggplot(CM_size_burntype, aes(as.numeric(trt_acres))) + geom_density(fill = 'grey90') +
  scale_x_continuous(n.breaks = 10)

# run the actual binomial model so you can look up values on the y axis
dat <- CM_size_burntype %>% 
  mutate(burntype = ifelse(any_broadcast_burn == T, 'Broadcast burn', 'Pile burn')) 

m1 <- glm(data = dat, formula = VIIRS_detects ~ log(trt_acres) * any_broadcast_burn, 
          family = 'binomial')
summary(m1)

## 
## Call:
## glm(formula = VIIRS_detects ~ log(trt_acres) * any_broadcast_burn, 
##     family = "binomial", data = dat)
## 
## Coefficients:
##                                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                            -3.3515     0.7014  -4.779 1.77e-06 ***
## log(trt_acres)                          0.5774     0.1356   4.259 2.05e-05 ***
## any_broadcast_burnTRUE                  1.0729     1.0735   0.999    0.318    
## log(trt_acres):any_broadcast_burnTRUE  -0.1319     0.2106  -0.626    0.531    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 275.10  on 204  degrees of freedom
## Residual deviance: 238.72  on 201  degrees of freedom
## AIC: 246.72
## 
## Number of Fisher Scoring iterations: 4

predict(m1, 
        newdata = data.frame(trt_acres = c(25, 50), any_broadcast_burn = c(F, T)),
        type = 'response')

##         1         2 
## 0.1835056 0.3692242

m2 <- glm(data = dat, formula = VIIRS_detects ~ (trt_acres) * any_broadcast_burn, 
          family = 'binomial')
summary(m2)

## 
## Call:
## glm(formula = VIIRS_detects ~ (trt_acres) * any_broadcast_burn, 
##     family = "binomial", data = dat)
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -1.2512806  0.2645627  -4.730 2.25e-06 ***
## trt_acres                         0.0017479  0.0005323   3.284  0.00102 ** 
## any_broadcast_burnTRUE            1.0683878  0.3582930   2.982  0.00286 ** 
## trt_acres:any_broadcast_burnTRUE -0.0015726  0.0006651  -2.365  0.01805 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 275.10  on 204  degrees of freedom
## Residual deviance: 254.75  on 201  degrees of freedom
## AIC: 262.75
## 
## Number of Fisher Scoring iterations: 4

predict(m2, 
        newdata = data.frame(trt_acres = c(25, 50), any_broadcast_burn = c(F, T)),
        type = 'response')

##         1         2 
## 0.2301291 0.4565784